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ABSTRACT 


In this study, we present an algorithm for system identification for systolic arrav 
implementation. With this schema, discrete samples of input and output data of a sys- 
tem with uncertain characteristics are used to determine the parameters of its model. 
The identification algorithm is based on recursive least squares, QR decomposition, and 
block processing techniques with covariance resetting. The identification process is 
based on the use of Givens rotation. Additionally, we want to address the following 
problems: how the round-off error propagates in time and the implementation in closed 
loop adaptive control. We will compare the implementation of fixed point arithmetic 
with the implementation of floating point arithmetic. This is primarily a theoretical in- 
vestigation to be conducted with computer simulations where numerical results will be 
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I. INTRODUCTION 


A. BACKGROUND 

When we design control svstems We are faced with the problem of identifying the 
plant to be controlled. In particular, we trv to determine the parameters of a math- 
ematical model (i.e., a linear differential or difference equation) which best fits the 
input-output data of the given plant. 

In some cases, insight of a model mav be obtained from the laws of Physics, 
Chemistry, etc. In most cases, this may not be possible due to the complexity of the 
physical factors involved. In these instances, it may be possible to derive the values of 
the parameters by observing the nature of the system’s response under appropriate ex- 
perimental conditions. This procedure is called “parameter estimation”. 

The problem of adaptively controlling systems with uncertain characteristics de- 
pends on identification of the unknown system parameters. In these cases, the parame- 
ters of the controller are computed on the basis of the current estimate of the dynamics. 
Figure 1 shows a general structure of an adaptive control system. 

The purpose of parameter estimation is to best fit a proper model to the input- 
output data of the system under investigation. The main issues are the following: 


]. Select an appropriate class of models, and 


2. Select an appropriate estimation algorithm. 


For the particular case of estimating the parameters of linear models, we search 
Within the class of models with a given fixed order for the one which minimizes a pre- 
diction error criterion. 

Firstly, suppose we wish to select a model of a system based on input-output 


measurements, in the form 
3) = 9 (16 + v(0) (1.1) 


where \(t) and g7(t) are determined on the basis of the output and input signals and @ 
is an array of unknown parameters to be determinated. The term v(t) represents noise 
or other modeling errors. Further information about this equation will be given in the 
next chapter. This class of linear models is preferable because of its simplicity and the 


amount of theory developed to analyze them. 


Secondly, we need to select an appropriate estimation algorithm. Although a wide 
choice exists, the most effective in terms of speed of convergence and accuracy is the 
recursive least squares algorithm [Ref. 1]. 


Identification 





Figure 1. The Adaptive Control System. 


The origins of the least-squares method can be traced back to Gauss as in [Ref. 2]. 
In its recursive version it has been formulated by several authors. The major drawback 
of recursive least-squares identification is its cost in terms of complexity, which might 
make it unsuitable for real time identification of a large number of parameters, since the 
size of the matrices involved grow with the complexity of the the system to be estimated. 
Although available microprocessors are effective for low order systems and slow sampIl- 
ing rates, more complex problems require improved capabilities. 

In this thesis we address the problem of implementing recursive least-squares iden- 
ufication using parallel processing techniques and systolic arrays. Particularly, the pa- 
rameter estimates are determined by the least-squares solution of a redundant number 
of linear equations obtained from the measured data. The techniques of solutions for 


this class of equations are based on QR decomposition, discussed in the next chapter. 


to 


B. SYSTOLIC ARRAY 

The idea of systolic arrav was developed by Kung and associates [Ref. 3] at 
Carnegie-Mellon University, and many versions of systolic processors are being designed 
by universities and industrial organizations. This subsection reviews the basic principle 
of systolic architectures and explains why they should result in cost-effective high per- 
formance, special purpose systems for a wide range of potential applications. 

A systolic array consists of a set of interconnected cells, each capable of performing 
simple operations. Since simple, regular communication and control structures have 
substantial advantages over complicated ones in design and implementation, the cells in 
a systolic svstem are typically interconnected with the immediate neighbors. In the 
systolic array we are considering there are two kinds of cells (boundary cell, internal 
cell). Each cell in the array is provided with local memory of its own, and it is connected 
only to its nearest neighbors. The array is designed such that regular streams of data 
are clocked through it in a highly rhythmic fashion, much like the pumping action of the 
human heart; hence the name “systolic”. Information in a systolic system flows between 
cells in a pipelined fashion. and communication with the outside world occurs onlv at 
the boundary cells. 

Basic principle of a systolic array is illustrated in Figure 2. Suppose each processing 
element in Figure 2 operates with a clock period of 100 ns. The conventional memory 
and processor organization in Figure 2a has 5 million operation per second. With same 
clock rate. the systolic array will result in 30 million operation per second performance 
provided the processing elements operate in parallel on pipelined data. The gain in 
processing speed has been increased six times. Being able to use each input data item a 
number of times is just one of the many advantages of the systolic approach. Other 
advantages include modular expandability, simple and regular data and control flows, 
use of simple and uniform cells and fast response time. 

Previous authors have presented parameter estimation algorithms using systolic 
arravs. The general idea has been to solve a system of linear equations in two stages: 


1}. Triangularization of the matrix of coefficients, 


2. Solving by successive substition. 


As explained in Chapter III the previous algorithms used a triangular svstolic array 
to triangularize the matrix, and a linear systolic array configuration to solve for the pa- 
rameters. The linear section requires operations such as divisions which are hard to 


implement by simple processor operations. Because of this, a new algorithm was 


developed. In our implementation, a second identical triangular array has been used 
instead of the linear systolic array. It is characterized by the fact that only orthogonal 
operations are involved, making the algorithm numerically more stable and easily 
implementable by simple shift and add operations. Also this algorithm needs two dif- 
ferent cells (boundary and internal); previous algorithms needed four different cells (two 
for a triangular array, two for a linear array). Although more total cells are required in 
this implementation, the cost of additional cells in a VLSI scheme is considered to be 


nunimial. 





Figure 2. The Concept of Systolic Processor Array. 


In implementations based on fixed point arithmetics, vector rotations necessary for 
the QR factorization can be implemented by a CORDIC algorithm, based on simple 
shift, add operations. 

The use of fixed point versus floating point arithmetic is considered during this in- 
vestigation. Because fixed point operations are based on simple shift functions and finite 
registers, Which are simple to implement, it seems advantageous to use fixed point val- 
ues. However, since input and output data do not naturally appear as integer values, 


there is concern over loss of accuracy due to necessary scaling and truncation. 


This research report is divided as follows: Chapter II discusses the methods of svs- 
tem identification, i.e., solution of systems of linear equations, QR decomposition, and 
recursive least squares algorithm, block processing and covariance resetting. Chapter 
III discusses the Givens rotation, the CORDIC technique and implementation of the 
systolic arravs. Chapter IV presents the simulation results, and Chapter V shows the 
final conclusions. A listing of the computer program is used to simulate the systolic 


arravs is found in the Appendix B. 


Il. MODELS FOR SYSTEM IDENTIFICATION 


A. LINEAR SYSTEM MODELING 


Suppose we wish to identify a model of a system based on input-output measure- 


ments. As shown in Figure 3 consider a system with a single input u(t) and single output 


y(t). 





Figure 3. Simple Linear System. 


If we consider a linear difference equation to model the dynamics, we can write 
KO + aylt— D+... + avi — n) = but — I) +... + 8,,u(t — ni) + (0) (2a) 
or, equivalently 


yD) = -ayt — WI —.... — ay(t— an) + butt — 1) +... +b, te(t — m) + v(2) 


“_~ 
ty 
to 

— 


where v(t) shows noise or other modeling errors, a,'s and 5, ’s are real constants, and the 
equation is of nth order. Since the output depends not only on the input but on the 
previous output values, this discrete system is known as recursive system. Equation (2.2) 


can be written in the matrix form 


t=) 

ii) 
Y(t) = Lair, ans Depo bm] (2.3) 

u(t — 1) 

u(t — m) 


where @ represents the parameter vector 

07 = 2 
Ser ag GeO pres Opn l (2.4) 

and y(t) represents the regression vector. 


—y(t— 1) 


- —y(r —n) > 5 
g(t) = nue) (2.5) 


u(t — m) 


By the above definitions we can write Equation (2.1) as 
¥(1) = 97 (0 + v(2) (2.6) 


The problem now is to estimate the parameter vector 8 from measurements of the 
sequences y(t) and g(r). Normally, if the number of samples of u(t) and y(t) (i.e., the 
number of equations) equals the number of unknowns (m+n), we can solve exactly for 
0. However, the number of equations is usually greater than the number of unknowns, 
since we tend to collect a large amount of data. For this reason, this system of equations 
in general does not have a solution. Ideally, in the noiseless case (v(t)=0) a solution 
exists, regardless of the number of equations. However, since noise and numerical errors 
are present, we look for the least squares solution of the svstem of equations, 1.e., for the 
solution which minimizes the error. 

If we take \ samples of Equation (2.6) and do not consider noise or other modeling 


errors We can write the system of N equations in matrix form as 


¥(2) gp" (t) 
x(¢— 1) gy (t—1) 
Bs 9 (2.7) 


ye-A)Y} | eUG-N 


or. equivalently 


b= AO (2.8) 


where 4 € 2,6 € BY, be BY. When N=M and A is full rank and invertible, we can 


solve uniquely for 6 as 
6=A'b (2.9) 


However, in signal processing applications we often face the case of M> N (i.e., more 
equations than unknowns), and the solution of (2.9) is defined in the least square sense 


by minimization of the error 

(9) = {140 — Bir (2.10) 
Therefore, the least squares solution 0, of (2.7) is implicitely defined as 

46, — Bll’ = mingl|46 — Bi (2.11) 


The least squares solution always exists, althought it might not be unique. We can 
solve Equation (2.11) by pseudoinverse which however, involves matrix inversion. An- 
other method is to triangularize A in Equation (2.8) and solve for @ by successive back 
substitutions. 

When M > N we meet the dilemma of triangularizing an array. The solution is to 
triangularize the upper part of the A matrix, leaving one or more rows of zeros at the 
bottom as a result. This will allow us to solve for 0 in the least square sense as indicated 


in Equation (2.11). This is known as QR decomposition. 


B. SOLUTION OF THE LEAST-SQUARES PROBLEM USING QR 
DECOMPOSITION 

A numerically attractive method for determing the least squares solution of a sys- 
tem of linear equations is provided by the QR decomposition of a matrix [Ref. 4]. 
Consider again Equation (2.8) with M > N. It can be shown that we can always 


factorize A as 
ACh (2) 


where A is a M x N matrix, Q is an M x M orthogonal matrix such that Q7Q = J and 
R is defined as 


ee ant (2213) 


where R is an N x N upper triangular matrix and 0 represents null matrix. It is not re- 
strictive to assume the diagonal entries of R to be non negative. Now we write again 


Equation (2.9) 

Ab=b 

with the equality intended in least squares sense, and multiply both sides by Q7 
Q’QRO=O'b 

sinee.O7O = | 

R0= 7b 


and, therefore. 





R 
| = A (2.14) 
0 2 
It is now easv to see that the error e(@) is given by 
46 — BI? = RO — Byll? + NBII” (2.15) 
which is minimal when 
RO=8, (2.16) 


since f, is independent of 8. In particular, if A is full rank and M > N, then the solution 
is unique, R is invertible and we can compute the least squares solution of (2.7) and (2.8) 


as 


Now the solution of Equation (2.8) in the least square sense can be computed in 
the same manner as the solution of a system of N equations and N unknowns in 
Equation (2.16). 

There is a number of methods available that can be used to compute the upper 
triangular matrix R (e.g., Householder transformation, Gram-Schmidt procedure, Giv- 
ens rotation). We will concentrate on Givens rotation. The Givens rotation is an at- 
tractive method for systolic array implementation since it successively manipulates two 
adjacent rows of the matrix at a time. CORDIC techniques can be used to implement 
these rotations and are discussed in the next chapter. 


C. RECURSIVE LEAST SQUARES ALGORITHM 

It is possible to estimate the values of the model parameters @ by observing the 
nature of the system’s response under appropriorite experimental conditions. The idea 
for the recursive identification problem is to compare the output with that of an adjust- 
able model, and update the parameters until the error between the outputs of the model 


and the system is minimized as in Figure 4. 


oy Suen 
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Figure 4. The Minimization of Error in the System. 


As seen in the previous section, assumption of a linear time invariant model leads to the 


relation 


H(2) = 97(0 (2.17) 


with @ the parameters to be estimated, and (rt) a measurable sequence. A recursive least 
squares procedure can be determined by defining 0, the estimate of @ at time t, as the 


vector which mininuzes 


tT] 
(0) = > [s) — 9290 |? (2.18) 


Particularly, 8,1s the least squares solution of the equations 


o7(t-1) r= 1) 


N 


A(t — 1)8, = M(t — 1) (2.19) 


a> 
= 

ll 

| 


| 670) ¥(0) 
It is possible to show that 6, can be recursively computed as [Ref. 2] 


ING) — O7(2) 


ey = 0 P(r — 1) - . (2.20) 
1+ 9 (1)P(t— I)e(s) 
Where P(t) = (4*()A(a))"' satisfies the recursion 
Pit— Yeo(ne (P(t — 1 
ee). eee (2.21) 


1+ 9 (P(t Lex) 


So far. we have neglected the problem of existence of the solution. The matrix A(t) 
can be singular, and then P(t) might not exist. Furthermore, we can not initialize A(-1) 
= 0 because, in this case this would vield P(-1) undefined, and we would not know how 
to start the sequence P(t). A solution to this problem is to incorporate initial conditions 
into A(t) so that P(t) can be computed at each t. Therefore, we modify the definition 
of 6, in Equation (2.19), so that initial conditions of the matrix P(t) can be accounted. 
Let A( —1) be any arbitrary matrix (it must be full rank). For convenience let A( —1) 
be square. with det (A( —1)) # 0, and define P(-1)= (4"( —1)A(—1))-'. This can be done 


by defining a new error criterion 


t—1 
(8) = > |x(s) — 07 (=) |? + (8 — 6)" 14( —)A( -1)) "0 — 6) (2.22) 


ay = 


which takes the initial estimate @, into account. This new formulation leads to the 


matrix A(t) as 
7 (0) 
¢7(¢-1) 


(2.23) 


We choose A(-1) = o,/, where o, > 0 is some arbitrary constant, and | the identity 
matrix of appropriate dimensions. In this way the matrix P(t) is always defined since 
AT(t)A(2) is never singular. Then, by algebraic manipulations, Equation (2.19) can be 


Written as 


o7(t— 1) Hea 1) 


AER A( —1)85 


The solution to (2.24) is computed recursively using (2.17) and (2.18) and the appropri- 


ate initial conditions. 


D. BLOCK PROCESSING AND COVARIANCE RESETTING 
From the previous considerations, the algorithm was described that allows us to es- 
timate the parameter 6 recursively. By using Equation (2.20) and P(t) = (A(2)7A(d)7) we 


obtain 
Pit) = (ATA)! = (+ > oO) (2.25) 
i=0 


t 
In general, the term .g(/)e7(i is likely to grow with time, so that P(t) — 0 as t ~ 00 


{Ref. 2]. Therefore, the algorithm loses sensitivity as t increases, and later values of @ 


may not be as accurate as earlier values, especially if our model changes with time. 
There are two possible solutions to this problem: 
1. The use of a “forgetting factor’, and 


2. The covariance resetting approach. 


By the forgetting factor approach we minimize the error 


Tt 


e()I? = ace = $7(k)O 44)? + o8l|O 41 — ll” (2.26) 


k=0 


where 0</<1 is the forgetting factor. This has the effect of assigning a higher weight 
to more recent data. 


The covariance resetting approach, on the other hand, divides the time scale into 
segments of equal and fixed length N as in Figure 5. At the end of each time block, 


KN 


Solve for 6 





Figure 5. Covariance Resetting Approach. 


we reset the covariance matrix P(t). Although it can be reset at any time, for conven- 


ience we choose the end of cach block, and P(t) now becomes 


-1 : 

Oo 1 r=kKN-1 k=0,!,2,... . 

fos : . (2.27) 
Equation(2.21) otherwise 


In particular, at the beginning of each interval the systolic array is initialized as o,/, and 
at the end of the interval (at KN+N-1 for the kth interval) the data are exited from the 
array and the system is solved to obtain Oe This estimate is used as an initial con- 
dition for the next time block. 

In an adaptive control context it has been shown [Ref. 5] that an external input u(t), 
sufficiently rich in frequency (n sinusoids), together with blocks I, of sufficent length N 


provides a guarantee for a consistent estimation as 
6— 0" 


where 6* represents actual parameters. The effect of various lengths of N will be inves- 
tigated later. 

If we apply the considerations given to the general case, we can see that the pa- 
rameter estimates at the end of each time block are related by the equation 


O(A+ IN - 1) WA + DN=1) 
o"(kN) Qty = (kN) (2.28) 
cay! oes 


A 


The systolic array implementation is based on this equation. In particular, 9,_),y 
is computed from Equation (2.28) at the end of each time block by making the leftmost 
matrix upper triangular. 

Although, the algorithm presented in this Chapter assumes a single input. single 
output (SISO) plant. it can be extended to the multiple input-output (MIMO) plant 
[Ref. 6]. Additionally, we assume the plant to be causal and of known order. 
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Hl. SYSTOLIC ARRAY IMPLEMENTATION 


Now, the problem is to compute the solution of Equation (2.28) using systolic ar- 
rays. In particular, Equation (2.28) is appropriate for parallel implementation. As de- 
scribed in [Ref. 5], the identification problem can be constructed again into a set of linear 


equations as 
RGusiyy = By (3.1) 


where R, is in upper triangular form. We can compute 0,.)y by using two processors 
in cascade: one to compute R, and f,, , and the second to compute @ from Equation 
(3.1). Notice that Equation (3.1) does not require any explicit matrix inversion since 
R, is in upper triangular form. 

As seen in Equation (2.28), we initialize the array at the beginning of each time 
block such that 


Ro = Gol 


By = 6,9; 


This has the effect of initializing the R, matrix in (3.1) to an upper triangular form. 
Then, at each discrete time t, an array of data g(r) and y(t) will be passed to the systolic 
array. The task of the array is to re-triangularize the data at each clock pulse, so the 
matrix remains in the form prescribed by (3.1). It is then a simple matter to solve for 
6... - This technique is based on QR decomposition as the means to triangularize the 
data array. 

The value of a, relates to the confidence we have in the initial estimates of 0. AS 
seen in Equation (2.26) a larger value of a, , will increase the confidence on the initial 


estimate of 6. Equations (2.24-2.28) show the role that o, plavs. 


A. GIVENS ROTATION 

The orthogonal triangularization process may be carried out by using various 
techniques. One of the techniques is the Gram-Schmidt orthogonalization procedure 
that provides the mathematical basis for the pipelined lattice predictor. Another pow- 


erful technique of triangularization by QR decomposition is provided by the Givens 
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rotation. The reason for this choice is the fact that the Givens rotation operates on pairs 
of adjacent rows, making it suitable for systolic array implementation. The object is to 
combine two adjacent rows in the matrix, forcing zeros in the appropriate positions so 


to obtain an upper triangular data matrix. Figure 6 shows an example of 


triangularization by successive Givens rotations, 





Figure 6. Example of Row Operations. 


Basic idea: Consider any two rows of data as 


a aK2 G80 arn 


Upp Akpy2 cee htt in 


We can see these two rows as a sequence of vectors in Figure 7. 

If we rotate all these vectors by an appropriate angle « then a,,,, 0 since the 
leftmost vector becomes parallel to the horizontal axis. This is accomplished by the 
following operation 


oY Ay y 


COS & Sin @ ar 1 ar 2 opooan Qk ny any ax 2 snocoo Qk n 








Se NA Gre Vale ie eo Ahn 0 gyi 2 ve Ak 1,n 


The value of the angle « can be determined using simple trigonometry 


Ak) : Qe+t 1 
cos © = sin a = 


2 2 2 2 
VN U1 F %H41,1 V G1 F Ake) 


This same rotation « is then applied to the remaining (x,y) vectors in the affected rows 


(3.2) 


(1.e., rows k+ 1, k) to ensure consistency. Next, the sequence of rotations is repeated for 
the reniainmimg pairs of rows (i.e., rows kK mek n-l? kK), ken Ro 


order to force zeros in the correct locations that leave an upper triangular matrix. 
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Figure 7. Example of Vector Operations. 


We remember that the data matrix is initialized to an upper triangular form (,! ). 
Then, as we take samples from the signals in the plant, new values of @ (t) and y(t) are 
added to the matrix. By a sequence of Givens rotations, we can transform it into an 
upper triangular matrix. An example set of rotations is shown in Figure 8, where the 
operations performed by the systolic array at each clock pulse are illustrated. This is 
repeated until t = KN (i.e., at the end of the time block), at which time the parameters 
6 are solved for, the matrix is reinitialized, and the process is repeated. 


Basis of the Givens rotation is the matrix 


im 0 0 
Ap.qy=|90 rpg) 0 (3.3) 
Oy “6, 


associated to each pair of indexes p, q € (1,n+1), with J, and J, identity matrices of di- 
mensions (q-2) x (q-2) and (n+ 1-q) x (n+ 1-q) respectively, and r is a 2 x 2 matrix of 


the form 


c(p.9q) s(p,q) 


a —s(p,q)  c(p.q) 


(3.4) 
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The matrix Q(p,q) is an orthogonal matrix and has the property that it affects only two 


rows at a time, ie., rows q and q-l. Application of the transformation Q(p,q) to any 


matrix A of appropriate dimensions implies 


x 


a 


x 


x 


Q(p.q)A = | 
aX 


x 


x“ ON 


as 


where x indicates other elements of the matrix and z indicates 


fs 2 
Z=V/%G-1p t Gp 


ot) p(t) xt) 
ee? oamcaie 


a 


Figure 8. Example of Givens Rotation. 


In Equation (3.4) c(p,q) and s(p,q) are provided such that 


C(P.9)Qg_y p + (P19) agp = Z 


C(P.9) Aq p — s(p.q)a, => 0 


(3.5) 





(3.6) 


In Figure 8, an example of application of the Givens rotation Q(3,4)Q(2,3)Q(1,2) to the 
left matrix results in the rotated matrix shown on the right (Ref. 7]. As seen before the 


Givens rotation requires addition, substraction, multiplication, division, and squares. 


Figure 9 illustrates the operations of each of the cells. Next, we will see how to perform 


the rotations, by using add and shift operations only in the CORDIC technique algo- 


rithm. 


Figure 9. 


EDGE EEEL 


(c out »S out ) 
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oa u/(sgrt(u +a )) 


S out = a/(sqrt(u7+ a”) 


a = sart(u*+ 6) 


endif 


INTERNAL CELL 





Givens Rotation Algorithm. 


B. THE CORDIC TECHNIQUE 

Here is a simple, accurate and relatively fast method for generating trigonometric 
functions. The technique proposed is based on the ingenious coordinate rotation digital 
computer (CORDIC) technique developed by J. E. Volder [Ref. 8]. 

The principle involved in CORDIC ts to rotate a vector, represented by its com- 
ponents X and Y, through a set of successive elementary rotations as seen in Figure 10. 


Each single rotation of the vector (X,Y) is computed at each step by a combination of 


simple add, subtract, and shift operations. 





Figure 10. Example of CORDIC Rotation. 


Consider a vector R represented by its components X and Y and rotate it through 
an angle + @. The results are the rotated coordinates \” and Y’ , which are related to X 


and Y by the following equations 


X'=Xcos@+¥siné@ 


(3.7) 
Y= Ycos? 2 isin? 


where the top sign refers to clockwise rotation. 


The CORDIC principle is based on performing a vector rotation in a sequence of 


angular steps , such that the sum of them equals to @, that is 


i) So teay am ee er amie (3.8) 
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If we define y = +1, then we can express 6 as 
n 

0= >) ym (3.9) 
t=0 


The problem is to choose the value of a, such that the equations above can be imple- 
mented with simple add and shift operations. This is possible if the values are chosen 
such that 


a,=tan (2) §=0,1,2...0 (3.10) 


Dividing both sides of Equation (3.7) by cos 0, gives 








OH I : 
COs 8, = AX;,+ Y; tan 6;= X45, 
aol 
Y", -* (3.11) 


The factors X”/ cos @, and ¥"/ cos @, are the rotated components of the vector (X,Y). The 
new vector is not only rotated, but also scaled by a factor 1/ cos 6,. If 6 is now replaced 
by o, = tan-'(2-') the new equations 
Nay =X; 4 V2 =X) 4 V7 = KN 

; (3.12) 
Peele = ey kT any 


The terms X, and Y’ indicate the components of the initial vector R. The term K;, is the 
factor 1/ cos @, by which X,., and Y,,, are larger than the components \” and Y’ of a 
vector which is only rotated. 

Rotations by angles smaller than 90°, the index (i) starts with 0,1,2,.... and contin- 
ues ton. The resulting values for 2~‘, a; and K, are given in Table 1. During each ro- 
tation the magnitude of R increases by K,=1/cose,. After n rotation steps the value 


of K, becomes 
l=n 


; 3 l ! | l 
K=| [k= Coe cone,  ‘cOson. cosa. (3.13) 
i=0 


| 


The value of K, thus increases with each rotation step and approaches 1.6468 for high 


values of i. 


Table 1. THE BINARY CORDIC CONSTANTS. 


9 ] 
1.0000000000 
0.5000000000 
0.2500000000 
0.1250000000 
0.0625000000 
0.03 12500000 
0.0156250000 
0.0078 125000 
0.0039062500 
0.001953 1250 
0.0009765625 


—, 


0 
1 
2 
S 
4 
5 
6 
7 
8 
9 
0 


—= 


x; (in degrees) 


45.00000000 
26.56505124 
14.03624340 
7.12501632 
3.57633432 
1.78991064 
0.89517384 
0.44761428 
0.22381056 
0.11190564 
0.05595300 


Kj 


1.414213500 
Se isaaZo 
1.629800596 
1.642484060 
1.645688908 
1.646492240 
1.646639215 
1.646743467 
1.646756030 
1.646759170 
1.646759954 





Since we assume the number of rotations to be constant, K, is also constant. And 
to correct for the encrease of the magnitude of R, the resulting VY, and Y, values must 
be divided by K, after the rotation operation is completed. Figure 11 illustrates a series 
of rotations for an arbitrary vector which we desire to rotate by 40°. Notice that in the 
figure, the desired angular value is approached quickly, because a, decreases by approx- 
imately half during cach step. The y sequence in the figure would be ( +1, -1, +1, 
+1, l,l, lel 

In an ideal case such as the one presented, the value of Y decreases during each 
step. However, this may not always be the case. For example, consider the vector (X,Y) 
= (5,2) anda = 12.5°. The first rotation in the CORDIC algorithm would be —45°, this 
gives us a rotated vector (X,Y) = (4.94,-2.12) and a» = —32.5°. Notice that the value 
of | }’| has actually increased. To make the algorithm more efficent, we modify the se- 


quence y to include the value zero. Whenever a rotation causes { }’] to increase, we do 


to 
tw 





Figure 11. Set of CORDIC Rotations. 


not perform.it, and we set the corresponding y, to zero. Then, we continue on with the 
next rotation value and repeat the process. For the example just mentioned, the next 


value to be tried would be 26.6°.. A CORDIC rotation algorithm can be seen in Figure 
12. 


C. SYSTOLIC ARRAYS 

In this section we will examine the parallel structure that will be used to solve the 
least squares algorithm described above. In particular we aim at a structure that accepts 
a sequence of regression vectors y() and signal y(n) as input and then outputs the es- 
timate for the parameter 0. Specifically, we are interested in a high performance parallel 
structure that can be implemented directly as a hardware device in order to deliver 
maximum throughput. Systolic arrays represent a structure suitable for these charac- 


teristics. : 
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* 


Initialize: If x>0O then x =x, y =y 
else X =-xX,Yy =-y 
if y> O then by= +1 
else py S| 


i=0 

while eee do 

if yi- 0) 2 xj > y, then 

Rit ST 

Yior = Yi 

G4, %= 0 

else e 

Kiet T= xX Ut Bi 2Yi 

Yivt = Yi- Bj2 x; 

if y >O then %,, 
else i,j 

endif 

endwhile. 





Figure 12. CORDIC Rotation Algorithm. 


A systolic array has simple and regular communication and control structures, and 
this is a substantial advantage over other designs and implementations. Additionally, 
we want to allow the computations to proceed concurrently with the input, tn order to 
maximize the throughput. This is known as pipelining [Ref. 3}. 

Figure 13 shows a typical system. As said before, the array is simply a network of 
processors that are regularly connected. The data is continuously “pumped” through 
this structure, thereby minimizing overall execution time, since all the processors work 
in parallel. 

It was shown that two processors would be necessary to solve Equation (3.1). 
Previously used configurations have consisted of a triangular array (as in Figure 13) to 


compute the upper triangular matrix, and a linear array to solve the system of equations. 
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Figure 13. Systolic Array. 


The entire systolic array is controlled by a single clock. Figure 14 shows the typical de- 
sign for the case d (dimension) = 3, where d shows the number of columns which enter 
into the system. In this thesis we will use the alternative configuration in which the 
linear section is replaced by a second triangular section identical to the first one. This 
new design will be discussed later in this chapter. Both designs use a single clock signal 
to control operations. Before continuing on to discuss the alternative design, we will 


review the structures of the triangular and linear sections. 


25 


Triangular Systolic 


Array 
Linear Systolic 


Array 





Figure 14. Systolic Array Procedure for Old Design. 


1. Triangular Array 
The triangular systolic array performs a sequence of Givens rotations. The 
CORDIC technique is a possible implementation for the rotations. The processors work 
simultaneously at each clock pulse. The data regression vector p(t) and output signal 
¥(t) are inputs to the top of the array, and rotations are calculated at each clock cycle. 

The triangular array consists of two types of cells: edge cells (or boundary cells) 
and internal cells. The edge cells are represented by the circles and the boundary cells 
are represented by the squares as seen in the Figure 13 and Figure 14. Figure 9 defines 
the operations of these cells. The edge cell computes the rotation parameters c and s 
as seen in Equation (3.2). Each cell of the triangular array stores an element of the up- 
per triangular matrix R(n) from Equation (3.1), and it is initialized to zero for internal 
cells and to o,/ for the edge cells as mentioned before. The function of each row of 
processing cells in the triangular systolic array section is to combine one row of the 
stored triangular matrix with a vector of data received from the above cells in such a 
way that the leading element of the received data vector is annihilated. The data vector 
so obtained is then passed downward on to the next row of cells. The boundary cells in 
each row of the section computes the rotation parameters and then passes them to the 
right on the next clock cycle. The internal cells apply the same rotation parameter to 
all other elements in the received data vector. 

A delay of one clock cycle per cell is incurred when passing the rotation pa- 
rameters along a row. That is why it is necessary to “skew” the input data as seen in 
Figure 14, so that the input data interacts properly with the previously stored triangular 
matrix. Because the cells are operating simultaneously. the data in the system at any 
tume t consists of values from (2n) different matrices. Figure 15 demonstrates this for a 
system with n=3. In this figure, we can see that at time (t+ 5), there are also values 
present from the five previous matrices (i.e., t+ 4,t+3.,......,t). In order to get all the cells 
in the array to a similar time state, the array would have to be clocked an additional 2n-1 
(five) cycles, feeding zeros as input where necessary. At the completion, all cells will be 
eucne same time (t+ 5) [Ref.. 7]. 

Note that at the same time the triangularization process is being carried out, 
the column vector f,, is also being computed by the rightmost column of internal cells 


using ¥(n) as its input. 
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Figure 15. Data Flow in Triangular Section. 


At the end of the triangularization period (N), we are ready to solve for @4.,,y. 
The data in this triangular array is clocked out to the next array section that will com- 
pute the parameters. 

2. Linear Array 

The linear systolic array has been used in previous implementations to solve for 
the estimated parameters. The linear section consists of one boundary cell and (d-1) 
internal cells as seen in Figure 14. 

The operation of the cells as they compute the parameters are shown in Figure 
16. Note that the cells are different from those in the triangular array, increasing to four 
the number of unique cells necessary in the combined system. 

It is shown in [Ref. 3] that the time required to solve for 0 using the linear array 
is equal to 2d. At the end of the period, the parameters are used as initial values for the 
triangular array, and the triangular section again begins the operation. 


We now replace the linear section with a second triangular section, and discuss 


the differences between these two design approaches. 





Figure 16. Definition of Cell Operation for Linear Section. 


3. The Use of a Second Triangular Array as the Solution Section 
An alternative implementation can be obtained by solving 0 as shown in 
Figure 17. 


Triangular 
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Triangular 
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Figure 17. Systolic Array Procedure for New Design. 
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In this implementation, at the end of the triangularization period, the data is 
passed from the first triangular section to the second triangular section in a reversed 
position. The second triangular array performs the same type of operations as the first, 
therefore, the cells are identical. 

The reason of using another triangular section is that, by proper combinations 
of rows, we can force Zeros into all elements of a given row but one, so that we can solve 
for each of the parameters. Also, the fact that orthogonal operations are used makes it 
more robust in the presence of numerical errors [Ref. 7]. To see how this works, consider 
an arbitrary set of equations in upper triangular form as in Equation (3.14). Now x, 1s 
simply solved as b,/a,. To solve for x,, we can force a,,; to zero by a linear combination 
of rows two and three. Then x, is found to be 3,/a,,. Similarly, by row operations on 
row 1,2 and 3 we can make a,, = a,, = 0 in row one, and x, is found to be 4,/a,,. This type 


of operations are exactly what the triangular array was designed to do. 


QQ, G2 83) 1% by 
O a3. 43] |x2]} = 1 Dy (3.14) 
0 0 Q33 X3 b; 


Figure 1§ shows these operations in matrix format, and as mentioned before 
data is in reversed position. Notice that the array is initialized to all zeros here, whereas 
the first triangular section is initialized to 6,4,, and ol. 

To understand how the system operates, recall the first triangular array at 
ume N. Now we must feed the data down into the second triangular section in an 
appropriate manner so that we can solve for the parameters. The same delay (one clock 
cycle per cell) in propagation of data applies to this triangular section as it does in the 
first section. That is why we must carefully choose when to sample the array in order 
to get the correct values with which to calculate the solution. 

Figure 19 illustrates the data flow for a simple system where d = 3. The input 
data is skewed as it was in the triangular section. It can be seen that the values a,;, 
a, , and a,, are available at times N+1, N+4, and N+7 respectively. Similarly, the 
values of b;, b,, and b, are available at times N+ 4, N+6, and N+ 8. In general, for any 
size system n, the coefficients @, and outputs } are available as shown in Table 2. Note 
from the figure that the coefficients are “picked off’ from the edge cells at the appropri- 


ate times, while the outputs are found in the rightmost set of internal cells [Ref. 3]. 
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Figure 18. Solution of System of Equations. 


This new design operates slower than the linear svstem seen in [Ref. 7]. The 
time to solve for @ in this design is equal to the time until 6, appears, which is equal to 
(3n - 1). This compares to 2n (or 2d) in previous implementations: hence, n-1 more 
clock cycles are required. 

On the other hand, simplification is gained in the manufacturing process since 
the number of types of cells is reduced. The tradeoffs to be considered are simplicity 


(cost) versus speed. 
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Table 2. TINS OF AVAILABILTY OF DATA. 
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IV. SIMULATION STUDY 


1. Model Equation 
In general, a linear model can be expressed by either the difference equation 
relating the input and output sequence, or im a regression form as mentioned in the 
previous Chapters. In the example below. we will consider the problem of estimating the 
unknown parameters both with and without noise present. 
In our simulation program, we will use the first order model which has the 


discrete time transfer function such that 





b, 
He) = =— a, (4.1) 
which corresponds to the linear difference equation 
HO + ay(t — 1) = by u(r — 1) (4.2) 


Where u(t) and y(t) are the input and output sequences respectively. Equation (4.2) can 


be expressed as 


67 =[a,.b,] @ 


Ly 
— 


oy () = ir 1), ur — 1) (4.4) 
This corresponds to the regression model 

H() = 8" 9(0) (4.5) 
If the model has noise we obtain 

(0) = 8" p(s) + Ha) (4.6) 


For proper system identification the input sequence must contain a sufficient 


number of frequency components. In our simulating program we used the sine wave 


ep ee eee é 
) + sin( 10 ) + sin( 10 ) + sin( 


mL 


10 





671 ) 


u(t) = sin( 10 


(4.7) 
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For identification purposes the input must be sufficently rich in frequencies to 
excite all modes of the system. If we have fourth order system we should have at least 
four different sinusoidal components. 

In our simulation model we chose the values of the parameters as a, = 0.5 and 
b,=2. In the parameter estimation problem, these values will be assumed to be un- 
known. 

2. Noise Model 

The noise term v(t) is a sequence of independent random variables with zero 
mean values. The noise is added to both the incoming signal u(t) and measured output 
y(t). 

3. Choice of Initial Values 

For a recursive algorithm, we need to initialize the systolic array with an initial 
parameter estimate o,6(0) and a,/ in Equation (2.28). Here a, 1s related to the confidence 
of initial condition of 6(0). If some prior information about @(0) is available, it should 
be used for determining proper values of 6(0). If no prior information is available we 
choose 7) = On 

4. Choice of Block Length 

In our simulation program we have used four different block lengths both 
without noise and with noise in order to test the convergence rates in the different cases. 
We chose N = 3, 7, 10, 15 values for the block lengths, where N identifies the block 
sizes. 

5. Floating Point and Fixed Point Operations 

During simulation, the use of floating point versus fixed point arithmetic is 
considered. Fixed point arithmetic operations are performed using simple shift oper- 
ations and finite registers. Because of this, they are simpler to implement than floating 
point operations. On the other hand, since input and output data values do not na- 
turally appear as integer values, there is some concern over loss of accuracy. The sol- 
ution to the latter problem is to scale all the numbers so that they stay within the limits 


of the fixed registers. 
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A. NOISELESS MEASUREMENTS 
This is an ideal case. Figure 20 through 27 illustrate the results of these simulations. 
In these figures the estimated parameters ( 0, and 0, ) are plotted atong the vertical axis, 
the block number is plotted along the horizontal axis. In the following we will discuss 
the results for floating point and fixed point arithmetic. 
1. Results Using Floating Point Arithmetic 
Figures 20 through 23 illustrate the floating point results for block lengths of 
N = 3,7, 10 and 15 respectively. As seen from the figures the parameters exhibit the 
fastest rate of converge when N=7, Even though the estimated parameter a, converges 
very fast at N=15, the other estimated parameter b, converges slowly in this block 
length. As it can scen from the Table 3, the case of N=7 requires minimum number 
of clock cycles, therefore in this case we should choose a block length of seven. 
2. Results Using Fixed Point Arithmetic 
This situation has been simulated by adding the random noise to the measure- 
ments and the computations, so to account for round-off errors. 
Figure 24 through 27 illustrate the fixed point results for the same conditions. 
Notice that the parameters converge as fast as floating point operations. As scen from 
Table 4, for this case again, N=7 requires minimum number of clock cycles. This indi- 
cates that the degradation due to the implementation is not dramatic for the fixed point 
processors. When We consider the simplicity of fixed point processors, this is a distinct 


advantage. 


Table 3. THE RATE OF CONVERGE FOR FLOATING POINT(NO NOISE). 


BLOCK LENGTH(N) BLOCK NUMBER TOTAL 
AT CONVERGE CGC aCrCLES 
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Figure 20. Floating Point Operation for N= 3 (No Noise). 
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Figure 21. Floating Point Operation for N=7 (No Noise). 
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Figure 22. Floating Point Operation for N= 10 (No Noise). 
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Figure 23. Floating Point Operation for N= 15 (No Noise). 
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Figure 24. Fixed Point Operation for N=3 (No Noise). 
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Figure 25. Fixed Point Operation for N=7 (No Noise). 
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Figure 26. Fixed Point Operation for N= 10 (No Noise). 
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Figure 27. Fixed Point Operation for N= 15 (No Noise). 
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Table 4. THE RATE OF CONVERGE FOR FIXED POINT(NO NOISE). 


BLOCK LENGTH(N) BLOCK NUMBER TOTAL 
AT CONVERGE CLOEGYELES 





B. NOISY MEASUREMENTS 
As mentioned before the noise term v(t) is a sequence of independent random var- 
lables with zero mean. We used a variance of 0.5 for these random variables. As ex- 
pected in a real system, the noise is added to both incoming signal u(t) and output y(t). 
1. Results Using Floating Point Arithmetic 
Figures 28 through 31 show the results for the four different block lengths. 
As seen from the figures, for N=3 the parameters are most affected. For N=7 or 
N= 10, the parameters are less affected by the noise. When the block length is equal to 
7, the parameters have converged reasonably well within eight blocks (56 cycles). For 
N= 15, we see that the parameters are least affected by the noise. 
2. Results Using Fixed Point Arithmetic 
Figures 32 through 35 show the results for fixed point implementation. Again 
we see that they exhibit similar performance as in the floating point cases. Effects of 
block length are almost the same as described previously. 
To simulate fixed point behavior we added random disturbances in the com- 
putations within the cells. In particular, the factors c and s are affected by round-off 
errors which we can simulate in this fashion. A similar approach has been taken in [Ref. 


9} in the analysis of a systolic array implementation of the projection operator. 
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Figure 28. Floating Point Operation for N=3 (With Noise). 
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Figure 29. Floating Point Operation for N=7 (With Noise). 
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Figure 30. Floating Point Operation for N= 10 (With Noise). 
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Figure 31. Floating Point Operation for N= 15 (With Noise). 
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Figure 32. Fixed Point Operation for N=3 (With Noise). 
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F igure 33. Fixed Point Operation for N=7 (With Noise). 
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Figure 34. Fixed Point Operation for N= 10 (With Noise). 
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Figure 35. Fixed Point Operation for N= 15 (With Noise). 


V. CONCLUSIONS 


In the previous chapters we explained the estimation of the model using a parallel 
structure and recursive least squares algorithm. In this chapter we will discuss the 
tradeoffs, advantages and disadvantages of the systems that we have investigated. 

The recursive algorithm, based on recursive least squares with covariance resetting, 
is redesigned in order to make it suitable for implementation on a VLSI chip and also 
modify the systolic array in order to reduce computing time and complexity. The dif- 
ference in adaptive control is that the estimator has to operate recursively on subsequent 
blocks of data, so that the estimated parameters converge to the respective correct val- 
ues. This convergence is guaranteed by initializing each estimation with the parameter 
values from the previous one, and by persistency of excitation of the external inputs. 

We simulated several possible conditions to see the response of the parallel algo- 
rithm. Almost in each condition the parameters converged close to their own values, 
even in the presence of noise. 

This parallel algorithm is different from other similar algorithms in the fact that 
we replaced two different processors by two identical ones. There was a total of four 
different computing cells (two for triangular section, two for linear section) in the old 
design. For the new design we need only two different cells for the whole system as ex- 
plained in Chapter III]. The new design requires more total cells than the old design. 
For a fourth order system (d=4), new design requires 10 additional cells from the for- 
mula of _ dd +1). However, additional cells are not expensive in a VLSI schema. 

As described in Chapter III additional (d-1) clock cvcles are required to operate the 
second triangular section. That is why a second triangular section 1s a little bit slower 
than the linear section. 

In this thesis we also compared floating point operations with fixed point oper- 
ations. For floating point operations we saw that the convergence rate increases with 
the block length. Also when the block length is small, the identification is more sensetive 
to the presence of noise. 

We got almost the same results for fixed point operations. If we consider the 
simplicity of the fixed point processor, a significant advantage to use the fixed point 
arithmetic is due to its simplicity. 
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APPENDIX A. COMPUTER PROGRAM I 


A. PURPOSE OF THE PROGRAM 


Givens rotation algorithm. There are two classes of cells. Therefore, two different types 
of subroutines are used. As explained in Chapter III this upper triangular matrix again 


convert to another upper triangular matrix and meanwhile the unknown parameters are 


This program converts a given data matrix to an upper triangular matrix by using 


computed. 


matrix. If run this program, should be extended virtual storage capacitv to 1500K in 


All elements of the matrix are given interactively. This program can solve 49 x 50 


advance due to the large arrav. 


(G2) ($2) (3) 1) 2) (i a) 


wkiekioetees VARTABLE DECLARATION titerereteied 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 A,C,S,U,LU,AOUT,AA,ED, IN, AAA,CC,SS,UU,AAOUT, LLU,EED, 


Sur dC EIN We Be, IsS 


DIMENSION AC50,50) C(O: 50, 0:50) S(O: 50 70:50), UC 0350, 0. 50)R 


+LU(0: 50,0: 50) ,AOUT(0: 50,0: 50),4A(50,50),ED(0: 50,0: 50), 
+IN(0: 50,0: 50) ,AAA(50,50) ,CC(0: 50,0: 50) ,SS(0: 50,0: 50), 
+UU(0: 50, 0:50), LLU(0: 50,0: 50) ,AAOUT(0: 50,0: 50), 
+EED(0:50,0:50), IIN(0: 50,0: 50) , THETA(50) 


INTEGER 1,344, 


Tere seve Tere tee de Tere sete FeF8 Te Te Tee VE NTE TS ITER TERETE RTC TEN 


severe ewe VARTABLE DEFINIAT ION 8% eevee 






ACI,J) = DATA MATRIX 

ED(I,J) = FIRST TRIANGULAR MATRIX 
EED(I,J) = SECOND TRIANGULAR MATRIX 
THETA(K) = UNKNOWN PARAMETERS 


devevedederededesedetewede vote tevedeve seve weve te tose tee fete vee eee ese 
detededeteeetene eres T ATA ENTRANCE wvtrevedededs dedevedededs 
PRINT *,'ENTER THE NUMBER OF COLUMNS M=' 
READ(5,*)M 
PRINT *,'ENTER THE NUMBER OF ROWS N=' 
READ(5,*)N 
PRINT *,' ' 
DO 1 I=1,N 
DORZ a= ier 
WRITE(6,3)1,J 
FORMAT('ENTER A(',12,12,')') 
READ(5,*) A(I,J) 
CONTINUE 
CONTINUE 
PRINT *,'DATA MATRIX' 
DO 4 I=1,N 


56 


WRETE(6,*)CAC1,J),J=1,M) 


4 CONTINUE 
PRINT *," ' 
DO 5 I=1,N 
DO 6 J=1,M 
AACT,J)=AC(N-(I-1),J) 
6 CONTINUE 
5 CONTINUE 
Do? I=1,2,3 
woes J=1,M 
AA(N+I,J)=0. 0 
8 CONTINUE 
7 CONTINUE 
DO 9 I=1,2,3 
DO 10 J=1,M 
AOUT(I,J)=0. 0 
10 CONTINUE 
9 CONTINUE 
DO 11 I=1,N+1 
Danie J) 
UCI, J)=AACI,J) 
12 CONTINUE 
ala CONTINUE 
(6 DIT TTT EIR ER ERITRERERETEERRIRER 
C deiaekekAET RST TRIANGULARIZATION** 27% 
DOmicrk—1N 


DO 14 J = KM 

DO 15 I = 1,N-K+2 

IF (J. EQ. K)THEN 

CALL EDGE(U(I,J),AOUT(I,J) ,AOUT(I+1,J),C(I,J),S(I,J)) 

ED(I,J)=AOUT(I,J) 

ELSE 

CALL INTERNAL(UCI,J),AOUT(I,J),C(I,J-1),S(I,J-1) , AOUT(I+1,J), 
Teves) ,CCL,JjgsC1,J)) 

UC(I-1,J)=LU(I, J) 

IN(I,J)=AOUT(I,J) 


ENDIF 
15 CONTINUE 
14 CONTINUE 
13 CONTINUE 
DO 20 I=1,N 
DO 21 J=1,N+1-I 
ED(I,J)=0. 0 
21 CONTINUE 


20 CONTINUE 
DOs22 J=M,2,-1 
DO 23 K=N+3-J,N+1 
DO 24 I=K,N+1 
ED(I,J)=IN(I,J) 
26 CONTINUE 
23 CONTINUE 
22 CONTINUE 
PRINT *,'FIRST TRIANGULAR MATRIX’ 
DO 29 I=N+1,2,-1 
WRITE(6,*)(ED(I,J),J=1,M) 
29 CONTINUE 
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PRINT *,' ' 
Cc devevevededeseveseitee DATA ENTRANCE ‘esvedededetesetecededecs 
DG 25 1=10N 
DO 26 J=1,M-1 
AAA(1,J)=ED(1+1,M-J) 
26 CONTINUE 
25 CONTINUE 
DO 27 J=M,M+1,2 
DO 28 I=1,N 
AAA(1,J)=ED(1+1,M) 
28 CONTINUE 
724) CONTINUE 
DO 30 I=1,2,3 
DO 31 J=1,M 
AAA(N+I1,J)=0. 0 
31 CONTINUE 
30 CONTINUE 
DOPS Zilles 
DO 33 J=1,M 
AAOUT(1,J)=0.0 
Ss CONTINUE 
32 CONTINUE 
DO 34 I=1,N+1 


DO 35 J=1,M 
UU(I,J)=AAA(I,J) 
35 CONTINUE 
34 CONTINUE 
(G dedevodevetovedededevededede fevetededede cece deve te deve de tote Te cede seed eters 
G dedededededesers SECOND TRIANGULARIZAT ION ttre 
DO 36 K=1,N 


DO 37 J = K,M 
DO Soule aN-kae 
IF (J. EQ. K)THEN 
CALL EDGE(UUC(1,J),AAOUI( 1,3) ,AAOUR( 1+] Gel SsCr. J) 
EED(I,J)=AAOUT(I,J) 
ELSE 
CALL INTERNAL(UU(I,J) ,AAOUT(I,J),CC(I,J-1),SS(I,J-1), 
+AAOQUT( I1+1,J), LLUC1,J).GUCT.)mecscl. J) 
UUCI-1, J)=ELUC 1d 
IIN(1,J)=AAOUT(I, J) 


ENDIF 
38 CONTINUE 
37) CONTINUE 


THETA(K)=IIN(2,M) /EED(2,K) 
36 CONTINUE 


DO 39 I=1,N 
DO 40 J=1,M-I 
EED(1,J)=0. 0 
40 CONTINUE 


39 CONTINUE 
DO 41 J=M,2,-1 
DO 42 K=M+2-J,N+1 
DO 43 I=K,N+1 
EED(1,J)=IIN(I,J) 
43 CONTINUE 
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42 
41 


44 


46 


45 


GG Gl G) 


CONTINUE 
CONTINUE 
PRINT *,'SECOND TRIANGULAR MATRIX’ 
DO 44 I=Nt+1,2,-1 
WRITE(6,*)(EED(1I,J) ,J=1,M) 
CONTINUE 
PRINT *,' ' 
PRINT *,'UNKNOWN PARAMETERS’ 
WRITE(* ,46) 
MORMAT(’ ' ,11Xp'K' , 13%) °THETA(K)') 
We 45) K=N,1,-1 
WRITE(6,*) N-K+1,THETA(K) 
CONTINUE 
STOP 


END 
Fereretedededevere rere de sete de Terese He HF Te II HIN HK TEI IE TIN 


ale lonteatentoats atealoaloatealoalsale alone ale aloalsale afeale ala wloalsatonfoatoatcntoatontlentaatos!snteate ol, ss 
ere er ere ete ee ee ere eee eee ee ee ee ee eee ee ee eel 


dietetics SUBROUTINE GROUPS FORsi#tidk 
detvtokdedststeveeeee CELLS FUNCTIONS iititeecddeded 


atoalentselen'selsateatan's alealon'selenisalen'selsclentsn'sa'se’on'selsalan'salsslsctentsclsaectealsn'snlectontec 
VEER EN CT PCED ES CRED OLED ECCS GOTT GL OS EN CD SEER CREEL GN CEE TEGO ELEC SEEDER STEVEN ED 


destededetsteiek eee BOUNDARY CELL ttttrstededetidededereds 
SUBROUTINE EDGE (U,AIN,AOUT,C,S) 
IMPLICIT REAL*8 (A-H,O-Z) 

REAL*8 U,AIN,AOUT,C,S 

IF(U. EQ. 0. AND. AIN. EQ. 0) THEN 

C=1.0 

S=0. 0 

AOUT=U 

ESE 

C=U/( SQRT( U**2. ODO+AIN**2. ODO) ) 
S=AIN/(SQRT( U***2, ODO+AIN**2. ODO) ) 
AOUT=SQRT( U***2. ODO+AIN**2. ODO) 
ENDIF 

RETURN 

END 


dedededededededetedetete tevet ete ee Tete TEETER ETE ETE RR EES 


sevevededsdevevereveve® INTERNAL CELL tecerevevevedsvevededetere 


SUBROUTINE INTERNAL (PU,AIN,PC,PS,AOUT,LU,LC,LS) 
PP ELCIT REAL 8 (A-H,O-2) 

REAL*8 PS,PC,PU,S,C,U,AOUT,LU,AIN,LC,LS 

LU=( -PU*PS )+( AIN*PC) 

AOUT=( PU*PC )+(AIN*PS) 

LC=PC 

LS=PS 

RETURN 

END 


Se a se a eciacae war ee seasie ere cis ami cin sie sede Sie ig oe ae ie Sioa SESE a Se 


of 


deaederededkweaedededkedeskeds He seve dkeaesedke dedeskeskeakeseveske seve deslenk keskese sk 

Sedevlevevedeveakaedevevesesesesedese dens sracsesone sede seak obsess aesesk a's 

verte THE RESULTS OF THIS PROGRAM rede 
NTER THE NUMBER OF COLUMNS M= 


NTER THE NUMBER OF ROWS N= 


NMoOmMaAM|AI| 


ENTER A( 1 1) 
i 
ENTER A( 1 2) 
D 
ENTER A( 1 3) 
3 
ENTER A( 2 1) 
4 
ENTER A( 2 2) 
5 
ENTER A( 2 3) 
6 
DATA MATRIX 
1. 0000 2. 0000 3. 0000 
4.0000 5. 0000 6. 0000 
FIRST TRIANGULAR MATRIX 
4.1231 5.3357 6. 5484 
0. 0000 0. 7276 1.4552 
SECOND TRIANGULAR MATRIX 
5.3851 4.0852 6. 6850 
0. 0000 0.5570 -0. 5570 
UNKNOWN PARAMETERS 
K THETA(K) 
1 -1. 0000 
2 2. 0000 


ENTER THE NUMBER OF COLUMNS M= 


4 

ENTER THE NUMBER OF ROWS N= 
3 

ENTER A( 1 1) 
3 

ENTER A( 1 2) 
7 

ENTER A( 1 3) 
ii 

ENTER AC 1 4) 
20 

ENTER A( 2 1) 
5 

ENTER A( 2 2) 
9 
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ENTER A( 2 3) 
16 
ENTER A( 2 4) 
21 
ENTER A( 3 1) 
6 
ENTER A( 3 2) 
8 
ENTER A( 3 3) 
10 
ENTER A( 3 4) 
22 
DATA MATRIX 
3. 0000 7.0000 
5.0000 9. 0000 
6. 0000 8. 0000 
FIRST TRIANGULAR MATRIX 
8. 3666 13. 6256 
0. 0000 2. 8884 
0. 0000 0.0000 
SECOND TRIANGULAR MATRIX 
21-8403 13. 7818 
0. 0000 2.0151 
0. 0000 0. 0000 
UNKNOWN PARAMETERS 
K THETA(K) 
1 Sie 77 
2 5. 8888 
3 14444 


tate 
ice 
10. 


20: 


ake 


ENTER THE NUMBER OF COLUMNS M= 


ENTER THE NUMBER OF ROWS N= 


ENTER AC 1 
ENTER A( 1 
ENTER AC 1 
BIER AC 1 
ENTER AC 1 
ENDER AC 2 
ENTER@SC 2 


ENTER A( 2 


1) 
2) 
3) 
4) 
5) 
1) 
2) 


3) 


0000 
0000 
0000 


6774 


. 6670 
- 2345 


3s BILAL 
5 SelyS 


2269 
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20. 
Bde 
pee 


S35 
ao7g2 
> Ble 


S55. 
. 6038 
ee 


0000 
0000 
0000 
4982 


2270 


5305 


1812 


ENTER A( 2 4) 
ENTER A( 2 5) 
3 

ENTER AC 3 1) 


9 
ENTER A( 3 2) 
Pal 
ENTER A( 3 3) 
42 
ENTER A( 3 4) 
57 
ENTER AC 3 5) 
48 
ENTER A( 4 1) 
12 
ENTER A( 4 2) 
19 
ENTER A( 4 3) 
23 
ENTER AC 4 4) 
45 
ENTER A( 4 5) 
58 
DATA MATRIX 
6. 0000 9. 0000 
42.0000 14. 0000 
9.0000 21. 0000 
12. 0000 19. 0000 
FIRST TRIANGULAR MATRIX 
44.9999 23. 5333 
0. 0000 22.9168 
0. 0000 0. 0000 
0. 0000 0.0000 
SECOND TRIANGULAR MATRIX 
85. 9883 69. 3000 
0. 0000 30. 5859 
0. 0000 0. 0000 
0. 0000 0. 0000 
UNKNOWN PARAMETERS 
K THETA(K) 
a 0. 5067 
2 0.4447 
3 -1. 6290 
4 1.7985 


14. 
Bh 
42. 


69) 
26. 
14. 


Bec 
=O) 


. 0000 
. 0000 
. 0000 
. 0000 


Sooo 
sgool 
25607 
~ 4541 


- 6206 
. 4896 
. 9060 
ous 


31. 
58. 


Se 
8h 
-14, 


Ve 
335). 


APPENDIX B. COMPUTER PROGRAM II 


A. PURPOSE OF THE PROGRAM 
This program computes the estimated weight vector of least squares for first order, 


second order or third order model which has the discrete time transfer function such that 


bez 


ee 
Z+ az +a,Z+ a; 


where the all coefficients are given interactively. In this program, addition to the previ- 
ous program is used another two subroutines to supply the gaussian noise which corre- 
sponds to parameter v(t). As seen in the plots this program computes the unknown 
parameters by converging. In order to run this program should be extended virtual 
storage capacity to 2M. 


(G eves eek V ART ABLE DECLARAT ION te versed secede 

REAL C,S,U,LU,AOUT,ED,IN,X,Y,SIG,UIN, THETA,V,R,BM,AAA,CC,SS,UU, 
+LLU , AAOUT ,EED, IIN 

DIMENSION C(0:50,0:50),S(0: 50,0: 50) ,U(0: 50,0: 50) ,LU(0: 50,0: 50), 
+AOQUT( 0: 50,0: 50) ,ED(0: 50,0: 50), IN(0: 50,0: 50),X(50) ,UIN(0: 50), 
+AAA(50,50) ,CC(0: 50,0: 50),$S(0: 50,0: 50) ,UU(0: 50,0: 50), 
+LLU(0: 50,0: 50) ,AAOUT(0: 50,0: 50) ,EED(0: 50,0: 50), 
+IIN(0: 50,0: 50) ,Y(-3: 50) , THETA( 0: 50) 

TVTECER 1.J,K.M,N,L,T 


Cc devedesedetesee eV ARTABLE DEFINIAT ION etvevesevede 
c ED(I,J) = FIRST TRIANGULAR MATRIX 
C EED(I,J) = SECOND TRIANGULAR MATRIX 
C etal ae le Ee De ae ate oa aie ce eee oe eiecastentefiacis vie cieaie Cia trea ispecies) 
c 
500 PRINT *, ENTER THE SIZE OF MATRIX N BY N N=?! 
READ *, N 
PRINT *,' N= ‘',N 
PRINT *,' ' 
PRINT *,'ENTER THE ORDERS OF DEN. OF DIFF. EQ. **NO=?' 
READ *, NO 
PRINT *,' NO = ',NO 
PRINT *,' ' 
PRINT *,'HOW MANY BLOCKS DO YOU WANT TO ITERATE **NBN=?' 
READ *, NBN 
PRINT *,' NBN = ' ,NBN 
PRINT *,' ' 
PRINT * "ENTER THE VALUE OF SIGMA **SIG=?' 
READ *,SIG 


PRINT *,'SIG = ',SIG 
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PRENT eee 
IF (NOSEQ 1) THEN 
PRINT *, "FORMAT OF 1ST ORDER 


PRINT *,'ENTER THE COEFF. OF 
READ *,Al 
PRINT *,'Al = ',Al 


PRINT *,'ENTER THE COEFF. OF 
READ *,B1 

PRINT *,'B1 = ',B1 
WRITE(*, 46) 


ELSEIF (NO. EQ. 2) THEN 

PRINT *,'FORMAT OF 2ND ORDER 
++A2*Y(T-2)' 

PRINT *,'ENTER THE COEFF. OF 

READ *,Al 

PRUNGT 3 A oa col 


PRINT *,'ENTER THE COEFF. OF 
READ *,A2 
PRINT *,'A2 = ',A2 


PRINT *,'ENTER THE COEFF. OF 
READ *,B1 

PRINT *,'B1 = ',Bl 
WRITE(*,47) 


ELSEIF (NO.EQ.3) THEN 

PRINT *,'FORMAT OF 3RD ORDER 
++A2*Y(T-2)+A3*Y(T-3)' 

PRINT *,'ENTER THE COEFF. OF 

READ *,Al 

PRINT *,'Al = ',Al 


PRINT *,'ENTER THE COEFF. OF 
READ *,A2 
PRINT *,'A2 = ',A2 


PRINT *,'ENTER THE COEFF. OF 
READ *,A3 
PRINT *,'A3 = ',A3 


PRINT *,'ENTER THE COEFF. OF 
READ *,B1 

PRINT * Gie= '. 8] 
WRITE(*,48) 


ELSE 


DIFF. EQ. ** Y(T) = B1*U(T-1)+A1*Y(T-1) 
la 


isle 


DIFF. EQ, #7 YC) = Si*U( t=) +Al-yCr-)) 


Ae. 


A2 ?' 


a2) 


DIFF. EQ. ** ¥(T) = BI*UCT-1)+A1*Y(T-1) 


ALee 


Oe 


ING ay 


Bilas 


PRINT *,'ERROR, PLEASE TRY AGAIN' 


PRINT *,'CHOOSE NO=1,2 OR 3' 
GOTO 500 

ENDIF 

NB=15*NO 

2=1 


ara) 


PI=4*ATAN(Z) 
P=PI/10 
BM=0 
R=0. 0707 
Tel 
Dedeve dete Fees Tete Fe Fe Te Fess Me Fe He Te Fe Fee ese Me FN Te Mee He eI Te HC 
wesecsdrsedeoeses ree DATA ENTRANCE tt teotee trsksect dese sescce 
DO 1 LF=2,NO0+2 
DO 2 KF=1,NO0+2 
IF( (LF. GT. KF). AND. ((LF-KF). EQ. 1)) THEN 
UC LF,KF)=SIG 
ELSE 
UC LF ,KF)=0. 0 
ENDIF 
CONTINUE 
CONTINUE 
IF(NO. EQ. 1) THEN 
Y(0)=0. 0 
UIN(0)=0. 0 
Y(1)=A1*Y(0)+B1*UIN(0) 
ESE bE GNOSEO. 2 HEN 
YC-1)=0. 0 
Y(0}=0. 0 
UIN(C0)=0. 0 
Y(1)=A2*Y( -1)+A1*Y(0)+B1*UIN(0) 
bESE 
Y(-2)=0.0 
Y(€-1)=0. 0 
Y¥(0)=0. 0 
UIN(0)=0. 0 
Y(1)=A3*Y(-2)+A2*Y( -1)+A1*Y(0)+B1*UIN(0) 
ENDIF 
DO 3 IT=1,NBN 
DO & J=1,NB 
CALL GAUSS( IX,R,BM,V) 
LE GNO EO mo) HEN 
U(1,1)=Y(J-1) 
UCL 2)=UIN(J-1) 
UGS) —VC J) 
UIN(J)=SIN( P*J)+SIN( P*¥2*J)+SINCP*5*J)+SIN( P*6%*J) 
Y¥(J+1)=A1*Y(J)+B1*UINCJ)+V 
BiSe (3 CNOVEO. 2) THEN 
UG = vGi-2)) 
He) -1) 
UC 1,3)=UINCJ-1) 
UC1,4)=¥CJ) 
UIN(J)=SINCP*J)+S INC P*2*J)+S INC P*5*J)+SIN( P*6*T) 
Y¥(J+1)=A1*¥(J)+A2*Y(J-1)+B1*UIN(J)+V 
PEs 
WG) == 3) 
0(1,2)=Y(J-2) 
0¢1,3)="C 1-1) 
UC 1,4)=UINCJ-1) 
U(1,5)=¥CJ) 
UINC J)=SINCP*J)+S INC P*¥2*J)+SIN( P*5*J)+S INC P*6%*J ) 
Y( J+1)=A1*Y(J)+A2*Y( J-1)+A3*Y(J-2)+B1*UIN(J)+V 
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dal 
10 


14 


Wz 


16 


18 
17 


BA 


ENDIF 
110) 5S) WUE 253 
DO 6 JJ=1,N+1 
AOUTCII,JJ)=0. 0 
UCN+II, JJ)=0. 0 
CONTINUE 
CONTINUE 
diekeokeeke PT RST TRIANGULARIZATION* Xe 
KARR TERRE EREE RUBE ECR 
DO 7 K=1,N 
DO 8 J1 = K,N+1 
DO 9 I = 1,N-K+2 
TF (1 EOS Kk HEN 
CALL EDGE( U(1,J1) ,AOUT( 1,31), AQUTCI+1,31) CCi ie Sei) 
ED(1,J1)=A0UT(I,J1) 
ELSE 
CALL INTERNAL(U(I,J1),AOUT(I,J1),C(I,J1-1),S(I,J1-1), 


+AOUTCI+1,31),LUC1,31),6(1,51),5C1.J1)) 


U-1, 31) =1Uct sa 
IN(I,J1)=AOUT(1I,J1) 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 
DO 10 I=1,N 
DOM Pies) 2—eaN plea 
ED(1I,J2)=0. 0 
CONTINUE 
CONTINUE 
IDLO) A HSSWheil 5 2g Pal 
DO 13 K=N+3-J3, Ned 
DO 14 I=K,N+1 
ED(I,J3)=IN(1I,J3) 
CONTINUE 
CONTINUE 
CONTINUE 
PRINT 3, 0 
sedevedededesedekeeke eT ATA ENTRANCE Ser tdrvedevedecrd ceded 
wietereteeF OR SECOND TRIANGULAR ARRAY **2re% 
DO 15 I=1,N 
DO 16 J4=1,N 
AAACT,J4)=ED( I+1,N+1-J4) 
CONTINUE 
CONTINUE 
DO 17 J5=N+1,N+2,2 
DO 18 I=1,N 
AAA(I,J5)=ED( I+1,N+1) 
CONTINUE 
CONTINUE 
DOPZO RI =I 253 
DO 21 J6=1,Nt+1 
AAA(N+I, J6)=0. 0 
CONTINUE 
CONTINUE 
DOR 22a lees 
DO 23 J7=1,Nt+1 
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ZS 
22 


25 
24 


28 
27, 
2 


30 
2S, 


33 


31 


35 
34 


AAOUT(1,J7)=0. 0 
CONTINUE 
CONTINUE 
DO 24 I=1,N+1 
DO 25 J8=1,N+1 
UU( I, J8)=AAA(I, J8) 
CONTINUE 
CONTINUE 
ede sese te dese Tee A AAAI HE ENTE KIA TIN ARIK 
wate SECOND TRIANGULARI ZATION* * 3% 
DO 26 K=1,N 
DOWZ7 J9 = KANti 
DO. 28 1 = 1,N-Kt2 
IF (J9. EQ. K) THEN 
CALL EDGE(UU(I,J9) , AAOUT(I,J9) , AAOUT(I+1,J9) ,CC(1I,J9) ,SS(I,J9) 


+) 


EED(I,J9)=AAOUT(I,J9) 
ELSE 
CALL INTERNAL(UU(I,J9) ,AAOUT(I,J9),CC(I,J9-1),SS(1,J9-1), 


tenowncI+i 39), LLU( 1,39), CC(1,J9),SS8(1,J9)) 


WUC T-1,J9)=LLUCT, J9) 
ITIN(I,J9)=AAOUT(I,J9) 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 
PRINT *,' ' 
Sesede devete ded de Bese ve seve Fe Sede Pose Be Fe Fee VETS Be BEE ETE EIS OT TE 
DO 29 I=1,N 
DO 30 J2=1,N+1-I 
REVOLT, J2)=0.0 
CONTINUE 
CONTINUE 
DORs J3=N+1 25-1 
DO 32 K=N+3-J3,N+1 
DO 33 I=K,N+1 
EED(1,J3)=IN(1,J3) 
CONTINUE 
CONTINUE 
CONTINUE 
DO 34 MM=1,N+1 
DO $5 LL=2 ,N+1 
UU( LL,MM)=EED(N+3-LL, MM) 
CONTINUE 
CONTINUE 
IF(NO. EQ. 1) THEN 
IF(UU(3,3).EQ. 0. AND. UU(3,2). EQ. 0) THEN 
THETA( 2)=0. 0 
ELSE 
THETA( 2)=UU(3,3)/UU(3,2) 
THETA( 1)=(UU(2,3)-( THETA( 2)*UU(2,2)))/UU(2, 1) 
WRITE(6,*) THETA(1) ,THETA(2) 
ENDIF 
ELSEIF(NO. EQ. 2) THEN 
IF(UU(4,4). EQ. 0. AND. UU(4, 3). EQ. 0) THEN 
THETA(3)=0. 0 
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46 


48 


aaaAn 


ELSE 

THETA( 3)=UU(4,4)/UU(4, 3) 

THETA( 2)=(UU( 3,4) -(THETA(3)*UU(3,3)))/UU(3,2) 

THETA( 1)=(UU(2,4)-(THETA(3)*UU( 2 ,3)+THETA(2)*UU(2,2)))/UU(2,1) 
WRITE(6,%*) THETA(1),THETA(2),THETA(3) 

ENDIF 

ELSEIF(NO. EQ. 3) THEN 

IF(UU(5,5). EQ. 0. AND. UU(5,4). EQ. 0) THEN 

THETA(4)=0. 0 

ELSE 

THETA(4)=UU(5,5)/UU(5,4) 

THETA( 3)=(UU(4,5)-(THETA(4)*UU( 4,4) ))/UU(4,3) 

THETA( 2)=(UU(3,5)-(THETA(4)*UU( 3, 4)+THETA(3)*UU(3,3)))/UU(3,2) 
THETA( 1)=(UU(2,5)-(THETA(4)*UU( 2 ,4)+THETA( 3 )*UU(2,3)+THETA(2)* 


+UU(2 ,2)))/UU(2, 1) 


WRITE(6,*) THETA(1),THETA( 2) ,THETA(3) , THETA(4) 
ENDIF 

ENDIF 

CONTINUE 

IF(NO. EQ. 1) THEN 

Y(0)=Y(NB) 

UIN(0)=UIN(NB) 

Y¥(1)=Y(NB+1) 

ELSEIF(NO. EQ. 2)THEN 

Y¥(-1)=Y(NB-1) 

Y¥(0)=Y(NB) 

UIN(0)=UIN(NB) 

Y¥(1)=Y(NB+1) 

ELSE 

Y(-2)=Y(NB-2) 

Y(-1)=Y(NB-1) 

¥(0)=YC(NB) 

UIN(0)=UIN(NB) 

¥(1)=YC(NB+1) 

ENDIF 
CONTINUE 
FORMAT('' ,17X,'THETA1' ,17X, 'THETA2' ) 
FORMAT('’ , 17M, THETA’, 17%," THETA? 18s THETAS. 
FORMAT('' ,17X, ‘THETA1’ ,17X,’ THETA2 -18X,’ THETAS . 19M ern AG" ) 
STOP 
END 
tevededovekioksevodededciedesesescdeietedesRitedcdedcdedtcsvetevecteded te 


vededevetoe ee SUBROUTINE GROUPS FOR*vdevevedex 
ARAMA CR EE KCELLS FUNCTION eaastecedcdesedettseds 


dedededevodededevere ee BOUNDARY CELLS ieverevcdstesedeerk 
SUBROUTINE EDGE (U,AIN,AOUT,C,S) 

REAL U,AIN,AOUT,C,S 

IF(U. EQ. 0. AND. AIN. EQ. 0) THEN 

C=1.0 

S=0. 0 
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90 


AOUT=U 

ELSE 

CALL GAUSS(1X,R,BM,V) 
C=U/(SQRT(U**2, ODO+AIN**2, ODO) ) 
C=U/(SQRT(U**2,. ODO+AIN**2, ODO) )+V 
S=AIN/(SQRT(U**2. ODO+AIN**2. ODO) ) 
S=AIN/(SQRT(U**2. ODO+AIN**2. ODO) )+V 
AOUT=SQRT(U**2. ODO+AIN**2. ODO) 
ENDIF 

RETURN 


END 
RABAT RACARERIEKICK KKK RICK KRENEK RIECK I REI 


sieddede ease tek TNTERNAL CELL eek irkkekkkcik 


SUBROUTINE INTERNAL (PU,AIN,PC,PS,AOUT,LU,LC,LS) 
REAL PS,PC,PU,S,C,U,AOUT,LU,AIN,LC,LS 

LU=( -PU*PS)+( AIN*PC) 

AOUT=( PU*PC )+(AIN*PS) 

LC=PC 

LS=PS 

RETURN 

END 

devedevedesedevesede Teves a sede ys clove see ve sess seve We se sesele rele 


devise SUBROUTINE GROUPS FOR‘ txdetet 
dededeakdeseskeeeeKGAUSSTAN NOLS Eazedsaratstdedes ceased 
SUBROUTINE GAUSS (IX,S,4M,V) 

A=0. 0 

DO 90 I=1,12 

CALL RANDU( IX, IY,Y) 

IX=IY 

A=A+Y 

V=(A-6. 0)*S+AM 

RETURN 

END © 


SUBROUTINE RANDU (IX,1Y,YFL) 
IY=1X65539 

Tey) 5.6,6 
IY=IY+2147483647+1 

YFL=1Y 

YFL=YFL*. 4656613E-9 

RETURN 

END 


sevedetevedededetetesevedesedede seven seve sede deter Wesese reek KK 
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